Problems with LASSO for variable selection in Synthetic Control

Using the passage of a drug decriminilzation law in Oregon as a case study, I look at how covariate adjustment in synthetic control can produce wildly different results.

Should you rely on LASSO for variable selection in Synthetic Control?

A recent working paper by Pickett, Hill, and Cowan (2026) on the problems of using only the outcome variable from untreated units without coviariates in synthetic control has reminded me of a recent example of why controlling for covariates can be extremely important when conducting a synthetic control analysis. Analysis of the impact of drug decriminalization in Oregon in 2020 provides a key example of 1) how results can vary greatly depending on which covariates are controlled for; and 2) the potential pitfalls of using an automated variable selection method like LASSO to construct the weights for the synthetic control.

Analysis of the effect of Oregon Decriminalization

In 2020, Oregon passed a law decriminalizing drug possession for small quantities of drugs. Since then, several papers have examined the impact of decriminalization in Oregon. The first paper, using the classic synthetic control method, found that the law increased overdose deaths by 23% relative to the synthetic control (Spencer 2023). However, a recent paper by Zoorob et al. (2024) found that the increase in overdose deaths attributed to the impact of decriminalization disappeared when controlling for the introduction of fentanyl into the Oregon drug market. Their results, using the matrix completion method from the excellent fect package, suggest that decriminalization had little impact on overdose deaths.

Andrew Wheeler also has some great Python and R posts on using LASSO for synthetic control and conformal inference with Oregon decriminalization as a case study. In this post I am going to use the approach shown in Wheeler’s blog posts - i.e. LASSO with the R package glmnet - to recreate the findings from the paper by Zoorob et al. (2024). Data on overdose death rates and fentanyl seizures comes from the replication materials for the paper by Zoorob et al. (2024) on Harvard Dataverse.

Loading in the data

#load in the packages
library(tidyverse)
library(glmnet) # for LASSO
library(readr)
library(meboot) #for maximum entropy bootstrapping
library(gt) #for nice tables


#load in the data, drop Washington, create unique period id
dat <-  read_csv("data_final.csv")

#data cleaning
glmnet_data <- dat|> 
  #drop Washington, years > 2020
  dplyr::filter(Year <= 2022 & ST != "WA") |> 
  group_by(State) |> 
  #create unique period id (30 periods)
  mutate(rowid = seq(from = 1, to = 30, by =1)) |> 
  ungroup()

Plot of Oregon Overdose Deaths

At first glance, there seems to be a clear rise in overdose deaths in Oregon around the time that the decriminalization law took effect in early 2021. In the plot below, I display the overdose deaths per 100k population for Oregon from 2008 to 2022. Overdose deaths exhibit a clear upward trend in the 2020s, but with a clear pretrend of rising overdose deaths since the second quarter of 2019.

#plot opiod overdoses per 100k for Oregon
overdose_plot <- glmnet_data |> 
  filter(State == "Oregon") |> 
  ggplot(aes(x = Year_half)) +
  geom_line(aes(y = od_death_rate_half), color = "blue", group = 1)  +
  geom_vline(xintercept= "2021_1" ,lwd=1,colour="gray",linetype = 2) + 
  labs(title = "Oregon Overdose Deaths", 
       y = "Overdose Deaths per 100k",
       x = "Period") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))



overdose_plot

An alternative explanation, as pointed out by Zoorob et al. (2024), is that overdose deaths rose as fentanyl was introduced into the Oregon drug market. In the plot below, I display the percent of drugs recovered that contained fentanyl as submitted to The National Forensic Laboratory Information System (NFLIS). As the plot shows, the percentage of drugs recovered containing fentanyl started to increase around the first quarter of 2020 before skyrocketing from 10% to around 45% from the first quarter of 2021 to the second quarter of 2022. Such a dramatic increase provides a very plausible explanation for the increase in the overdose death rate in Oregon around the time of the passage of the decriminalization law.

Plot of Fentanyl introduction into Oregon drug market

#plot fentantyl introduction for Oregon
fentantyl_plot <- glmnet_data |> 
  filter(State == "Oregon") |> 
  ggplot(aes(x = Year_half)) +
  geom_line(aes(y = nflis_fentanyl_percent_total), color = "blue", group = 1)  +
  geom_vline(xintercept= "2021_1" ,lwd=1,colour="gray",linetype = 2) + 
  labs(title = "Oregon Percent Fentanyl of drugs recovered", 
       y = "Percent Fentanyl",
       x = "Period") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  


fentantyl_plot

Synthetic Control Analysis

To see the impact of controlling for the introduction of fentanyl into the Oregon drug market, I will run two analyses: 1) synthetic control with glmnet and the weights chosen with LASSO; 2) synthetic control with glmnet and LASSO but with the percentage of drugs recovered that were fentanyl included as a predictor in the analysis.

Importantly, in the first analysis I will also include the percentage of drugs recovered that were fentanyl in Oregon as a potential predictor, but allow LASSO to select among the potential predictors as weights rather than forcing the variable to be included as in the second analysis.

#pivoting the data wide
wide_glmnet_data <- glmnet_data |> 
  tidyr::pivot_wider(id_cols = rowid, names_from = State, values_from = od_death_rate_half) |>  
  dplyr::relocate(Oregon,.after = rowid) 

#pre outcome
pre_data <- wide_glmnet_data |> 
  filter(rowid < 27) 

#post outcome
post_data <- wide_glmnet_data |> 
  filter(rowid >= 27) 

#grabbing fentanyl percent total as a covariate
covariates <- glmnet_data |> 
  filter(State == "Oregon") |> 
  select(rowid, nflis_fentanyl_percent_total)

#pre covariate
cov_pre <- covariates |> 
  filter(rowid < 27) 

#post covariate
cov_post <-  covariates |> 
  filter(rowid >= 27) 

#add to the pre data
pre_data <- left_join(pre_data, cov_pre, by = "rowid") |> 
  as.matrix()

#add to the post data
post_data <- left_join(post_data, cov_post, by = "rowid") |> 
  as.matrix()

For the synthetic control, I’ll predict the pre treatment overdose deaths per 100k for Oregon using the overdose deaths rates for the other states minus Washington (who implemented a similar policy around the same time) and include the percentage of drugs recovered that were fentantyl in Oregon as an additional predictor.

In the second analysis, I will force the inclusion of the percentage of drugs recovered that were fentantyl in Oregon using the penalty factor argument , which allows you to control the amount of penalization for each predictor -i.e. setting the percentage of fentayl as 0 to force no penalization and all other predictors set to 1 to force penalization.

For both analyses I will limit the weights to be positive by setting the lower limits to 0 and the upper limits to 1 and include an intercept. Then I’ll predict the post treatment counterfactual using the weights selected from the pre-treatment period.

Code to run the LASSO regressions

#run the LASSO regressions for the pretreatment period

#unadjusted lasso
set.seed(07102026)
lasso <- glmnet(x = pre_data[,3:52], y = pre_data[,2],family="gaussian",
                    lower.limits=0, upper.limits =1, intercept=TRUE,standardize=FALSE,
                    alpha=1, nlambda = 1000)


cv.lasso <- cv.glmnet(x = pre_data[,3:52], y = pre_data[,2],family="gaussian",
                          lower.limits=0, upper.limits =1, intercept=TRUE,standardize=FALSE,
                          alpha=1, nlambda = 1000)



#set up the penalty factor to include covariates
#The fentanyl percent in Oregon is the predictor in column 50, 
#with 0 forcing no penalization (inclusion) and 1 for all others. 
penalty_factor <- c(rep(1, 49),0)


#add in the penalty factor for the adjusted analysis
set.seed(07102026)
lasso_adj <- glmnet(x = pre_data[,3:52], y = pre_data[,2],family="gaussian",
                     lower.limits=0, upper.limits =1, intercept=TRUE,standardize=FALSE,
                     alpha=1, nlambda = 1000,
                penalty.factor = penalty_factor)


cv.lasso_adj <- cv.glmnet(x = pre_data[,3:52], y = pre_data[,2],family="gaussian",
                      lower.limits=0, upper.limits =1, intercept=TRUE,standardize=FALSE,
                      alpha=1, nlambda = 1000, 
                      penalty.factor = penalty_factor)

Examining the synthetic control weights

#extracting the wesights

#Weights for unadjusted
weights <- coef(cv.lasso, s = "lambda.min") |> as.matrix() 
weights <- as.data.frame(weights) 
weights <- weights |> format(weights$lambda.min, scientific = FALSE) |> 
  arrange(desc(lambda.min)) |> 
  mutate(lambda.min = as.numeric(lambda.min)) |> 
  filter(lambda.min !=0) |> 
  tibble::rownames_to_column(var = "Variable") |> 
  rename("Weights" = 2)


#weights for the adjusted
weights2 <- coef(cv.lasso_adj, s = "lambda.min") |> as.matrix() 
weights2 <- as.data.frame(weights2) 
weights2 <- weights2 |> format(weights2$lambda.min, scientific = FALSE) |> 
  arrange(desc(lambda.min)) |> 
  mutate(lambda.min = as.numeric(lambda.min)) |> 
  filter(lambda.min !=0)|> 
  tibble::rownames_to_column(var = "Variable") |> 
  rename("Weights" = 2)


#make into nice tables with gt
weights |>
  gt::gt() |>
  gt::tab_header("Weights Unadjusted") |> 
  opt_row_striping(row_striping = FALSE)

Weights Unadjusted
Variable	Weights
(Intercept)	2.50253727
California	0.27811205
Nebraska	0.11640906
New Mexico	0.10277543
Nevada	0.03522820
South Carolina	0.02816837

weights2 |>
  gt::gt() |>
  gt::tab_header("Weights Adjusted") |> 
  opt_row_striping(row_striping = FALSE)

Weights Adjusted
Variable	Weights
(Intercept)	5.34860990
nflis_fentanyl_percent_total	0.42641957
New Mexico	0.09885041

In the unadjusted LASSO synthetic control analysis, the weights selected at the lambda that minimizes Mean Squared Error (lambda min) include an intercept along with a combination of overdose deaths per 100k population from the states of California, Nebraska, New Mexico, Nevada, and South Carolina. Critically, even though the percentage of drugs recovered that were fentanyl in Oregon was included as a predictor, LASSO did not select this variable for inclusion in the analysis. This suggests that analysts need to be very careful when using LASSO for variable selection in synthetic control as LASSO may not select a variable that has a clear causal relationship in predicting the outcome.

Looking at the weights for the adjusted LASSO synthetic control analysis - i.e. forcing the inclusion of the percentage of drugs recovered that were fentanyl in Oregon - we can see that they include the percentage of fentanyl, an intercept, and overdose deaths per 100k population from the state of New Mexico.

Code to predict the synthetic counterfactuals

#grab lambda min
lambda_min <- cv.lasso$lambda.min
lambda_min_adj <- cv.lasso_adj$lambda.min

#Get our pretreatment period predictions using lambda min from cross validation or sub in for reproducibility
pre_pred <- predict(lasso, newx = pre_data[,3:52], s = lambda_min)

#Get our post treatment predictions using lambda min from cross validation 
post_pred <- predict(lasso, newx = post_data[,3:52], s = lambda_min)

#Get our pre treatment adjusted predictions use lambda min from cross validation
pre_pred_adj <- predict(lasso_adj, newx = pre_data[,3:52], s = lambda_min_adj)

#Get our post treatment predictions using lambda min from cross validation 
post_pred_adj <- predict(lasso_adj, newx = post_data[,3:52], s = lambda_min_adj)



qtrs <- glmnet_data |> 
  slice(1:30) |> 
  select(Year, Year_half)

pred_table <- rbind(pre_pred, post_pred) |> 
  as.data.frame() |> 
  rename(predictions = 1)

pred_table2 <- rbind(pre_pred_adj, post_pred_adj) |> 
  as.data.frame() |> 
  rename(predictions_adj = 1)


pred_table <- cbind(qtrs, wide_glmnet_data$Oregon, pred_table, pred_table2) |> 
  mutate(rowid = seq(from = 1, to = 30, by = 1)) |> 
  rename(Oregon =3) |> 
  relocate(rowid,.before = Oregon) |> 
  mutate(att = Oregon - predictions,
         att_adj = Oregon - predictions_adj)

Adjusted vs Unadjusted Prediction Plots

In the plots below, we can clearly see that adjusting for the percent of drug recoveries that were fentanyl in Oregon dramatically changes the predicted counterfactual. In the un-adjusted plot, the predicted overdose death rate is lower than the observed overdose death rate in the post treatment period, suggesting that the passage of the law decriminalizing drugs in Oregon led to higher oversdose deaths than what would have occurred if the law was not implemented.

However, when adjusting for the percent of drug recoveries that were fentanyl in Oregon, the plot reverses entirely. The predicted overdose death rate in Oregon is actually much higher than the observed death rate in the post treatment period. This suggests that overdose deaths in the post treatment period were actually lower than what would have been expected given the level of fentanyl in the Oregon drug market. For both plots, the fit between the predicted and observed overdose death rates appear similar in the pretreament period, indicating that pretreatment fit may not always be the best guide when conducting a synthetic control analysis.

#lets graph these bad boys

colors <- c("Oregon Observed" = "black", "Oregon Predicted" = "blue")

#Plot with ggplot
plot_unadjust <- pred_table |> 
  ggplot(aes(x = Year_half)) +
  geom_line(aes(y = Oregon, color = "Oregon Observed", group = 1)) + geom_point(aes(y = Oregon), color = "black") +
  geom_line(aes(y = predictions, color = "Oregon Predicted", group = 1)) + geom_point(aes(y = predictions), color = "blue") +
  scale_y_continuous(breaks = seq(from = 0, to = 30, by =5), limits = c(0, 30)) +
  geom_vline(xintercept= "2021_1" ,lwd=1,colour="gray",
             linetype = 2) + theme_minimal()  +
  labs(title = "Observed Oregon vs Synthetic Oregon (Un-adjusted)", 
       y = "Drug Overdose Deaths per 100k",
       x = "Period",
       color = "Legend") +
  scale_color_manual(values = colors) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

plot_unadjust

plot_adjust <- pred_table |> 
  ggplot(aes(x = Year_half)) +
  geom_line(aes(y = Oregon, color = "Oregon Observed", group = 1)) + geom_point(aes(y = Oregon), color = "black") +
  geom_line(aes(y = predictions_adj, color = "Oregon Predicted", group = 1)) + geom_point(aes(y = predictions_adj), color = "blue") +
  scale_y_continuous(breaks = seq(from = 0, to = 30, by =5), limits = c(0, 30)) +
  geom_vline(xintercept= "2021_1" ,lwd=1,colour="gray",
             linetype = 2) + theme_minimal()  +
  labs(title = "Observed Oregon vs Synthetic Oregon (Adjusted)", 
       y = "Drug Overdose Deaths per 100k",
       x = "Period",
       color = "Legend") +
  scale_color_manual(values = colors) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))


plot_adjust

While the plots above display much larger effects relative to the analysis done by Zoorob et al. (2024) with matrix completion, in both cases they show that adjusting for the introduction of fentanyl into the Oregon drug market reverses the treatment effect. This can clearly be seen from the plots in Figure 3 from Zoorob et al. (2024) below.

Figure 3 (from Zoorob et al. 2024)

Constructing Confidence Intervals for the ATTs

For good measure, I plot the ATTs - i.e. the observed overdose death rate for Oregon minus the predicted overdose death rate - from both the adjusted and unadjusted synthetic control models with confidence intervals constructed via maximum entropy bootstrapping using the approach outlined in Wong et al. (2023) for inference with LASSO using glmnet. The confidence intervals are too narrow for my liking and the method is too laborious to be practical in many cases, but it provides a little known alternative to conformal prediction or placebo tests for inference.

Bootstrapping for inference is probably a bad idea with LASSO and glmnet, but I like the applied aspect of constructing 1000 bootstrap samples with a time series resampling method and constructing 1000 counterfactuals for inference by refitting LASSO on each one. However, I would be very skeptical of using this approach in practice as it tends to produce confidence intervals that are asymmetrical and too narrow in many cases. If I was using this for a paper or formal policy analysis, I would be much comfortable using conformal inference or the approaches for inference outlined by Cattaneo et al. (2025) and used in their excellent scpi R package.

#code to get confidence intervals for the ATT

#function for effect intervals
bootstrap_cum_effects_fixed <- function(boot_pre_list, 
                                        boot_post_list, 
                                        boot_pre_outcome_column,
                                        outcome_post, 
                                        outcome_post_column, 
                                        lambda_min, 
                                        conf_int,
                                        penalty = NULL) {
  
  #Set the outcome columns, start and end columns
  outcome_pre_column <- boot_pre_outcome_column 
  start_column <- boot_pre_outcome_column + 1
  end_column <- ncol(boot_pre_list[[1]])
  
  if(is.null(penalty)){
  #apply the lasso to each bootstrapped sample to predict the treated with the non treated Columns
  lasso_list <- purrr::map(boot_pre_list, function(df) (glmnet(x = df[,start_column:end_column], y = df[,outcome_pre_column],family="gaussian",
                                                               lower.limits=0, upper.limits =1, intercept=TRUE,standardize=FALSE,
                                                               alpha=1, nlambda = 100)))
  }
  
  
  if(!is.null(penalty)){
    #apply the lasso to each bootstrapped sample to predict the treated with the non treated Columns
    lasso_list <- purrr::map(boot_pre_list, function(df) (glmnet(x = df[,start_column:end_column], y = df[,outcome_pre_column],family="gaussian",
                                                                 lower.limits=0, upper.limits =1, intercept=TRUE,standardize=FALSE,
                                                                 penalty.factor = penalty,
                                                                 alpha=1, nlambda = 100)))
  }
  
  
  
  #make predictions using the model with the old lambda for the post treatment period for each sample
  pred_sim <- purrr::map2(lasso_list, boot_post_list, 
                          function(.lasso_list, .boot_post_list) predict(.lasso_list, newx = .boot_post_list, s =lambda_min))
  
  #Calculate the effects for each predicted series for the post period (oregon post - diff)
  eff_sim <- purrr::map(pred_sim, function(diff) {(outcome_post[,outcome_post_column] - diff[,1])})
  
  #get the cumulative effects - cumulative ssum Observed Oregon minus Predicted Oregon
  cum_eff_sim <- purrr::map_dfc(eff_sim, function(cumulative){(cumsum(cumulative))})
  
  #bind all the rows together
  eff_sim <- do.call(cbind, eff_sim)
  
  #calculate the cumulative effects intervals
  #set the names for the effects table
  table_names <- outcome_post[,1:(outcome_post_column - 1)]
  
  #set the alpha
  alpha  = (1 - conf_int)/2
  alpha_lower = 0 + alpha
  alpha_higher = 1 - alpha
  
  #calculate the effect intervals
  eff_low = c(apply(eff_sim, 1, quantile, probs = alpha_lower))
  
  eff_high = c(apply(eff_sim, 1, quantile, probs = alpha_higher))
  
  average_eff = c(apply(eff_sim, 1, mean))
  
  #cbind to make table
  eff_table <- cbind.data.frame(table_names, eff_low, eff_high, average_eff)
  
  #calculate the cumulative effects intervals
  cum_pred_low = c(apply(cum_eff_sim, 1, quantile, probs = alpha_lower))
  #Grab upper prediction/effects interval
  cum_pred_high = c(apply(cum_eff_sim, 1, quantile, probs = alpha_higher))
  #Get the average of the cumulative effects
  average_cum_effect = c(apply(cum_eff_sim, 1, mean))
  
  #cbind to make table
  cum_effect_table <- cbind.data.frame(table_names,cum_pred_low, cum_pred_high, average_cum_effect)
  
  #put the outcomes into a list
  effects_list <- list(eff_sim, cum_eff_sim, eff_table, cum_effect_table)
  
  #set the names
  names(effects_list) <- c("effects","cum_effects", 
                           "effects_intervals", "cum_effects_intervals")
  
  #return the list
  return(effects_list)
  
}

#grab data for bootstrapping
oregon_sim <- glmnet_data |> 
  dplyr::select(rowid, State, od_death_rate_half, nflis_fentanyl_percent_total)



#Grab 1000 bootstrap ensembles from the outcome column of the ODR_sim_matrix
set.seed(12112025)

#bootstrap of the outcome column for all 49 donor states
oregon_boot1 <-meboot::meboot(as.matrix(oregon_sim$od_death_rate_half), reps=1000, trim = list(trim=0.10, xmin=0, xmax = NULL))

#bootstrap of the percent fentanyl for Oregon
oregon_boot2 <-meboot::meboot(as.matrix(oregon_sim$nflis_fentanyl_percent_total), reps=1000, trim = list(trim=0.10, xmin=0, xmax = NULL))



#Make ensembles into a dataframe
oregon_boot1  <- oregon_boot1$ensemble |> 
  as.data.frame() 

oregon_boot2 <- oregon_boot2$ensemble |> 
  as.data.frame() 


#Add in the state and period columns and clean names
oregon_boot1 <- cbind(glmnet_data$rowid,glmnet_data$State, oregon_boot1) |> 
  dplyr::rename(rowid =1, state = 2)|> 
  dplyr::relocate(rowid,.before = 1)


#Add in the state and period columns and clean names
oregon_boot2 <- cbind(glmnet_data$rowid,glmnet_data$State, oregon_boot2) |> 
  dplyr::rename(rowid =1, state = 2)|> 
  dplyr::relocate(rowid,.before = 1) |> 
  filter(state == "Oregon")



#grab only our bootstrapped samples and convert to a list for the pre period
oregon_boot1_pre <- oregon_boot1 |> 
  filter(rowid < 27) 

oregon_boot2_pre <- oregon_boot2 |> 
  filter(rowid < 27)

#make into a list for the pretreatment period
oregon1_boot_pre_list <- as.list(oregon_boot1_pre[,3:ncol(oregon_boot1_pre)])

oregon2_boot_pre_list  <-  as.list(oregon_boot2_pre[,3:ncol(oregon_boot2_pre)])

#add back in the column names to each object in the list and rename for pivoting
#for our entire ensemble series
oregon1_boot_pre_list  <- purrr::map(oregon1_boot_pre_list, function(bind) (cbind.data.frame(oregon_boot1_pre$rowid, oregon_boot1_pre$state, bind))) |> 
  #rename the columns
  purrr::map(dplyr::rename, rowid = 1, state = 2) |> 
  purrr::map(~tidyr::pivot_wider(.x, names_from = "state", 
                                 values_from = "bind"))|> 
  #reorder to make it easier for glmnet - Oregon in the second column
  purrr::map(relocate, Oregon, .after = rowid) |> 
  #convert to a matrix for LASSO
  purrr::map(as.matrix)


#Do the second covariate
oregon2_boot_pre_list <- purrr::map(oregon2_boot_pre_list, function(bind) (cbind.data.frame(oregon_boot2_pre$rowid, bind))) |> 
  #rename the columns
  purrr::map(dplyr::rename, rowid = 1,fentanyl = 2)|> 
  #drop the the rowid
  purrr::map(dplyr::select, -rowid)|> 
  #convert to a matrix for LASSO
  purrr::map(as.matrix)

#bind both lists together
oregon_boot_pre_list <- purrr::map2(oregon1_boot_pre_list, oregon2_boot_pre_list, cbind)



#do the post period
#Make it into a list
oregon_boot1_post <- oregon_boot1 |> 
  filter(rowid >= 27) 

oregon_boot2_post <- oregon_boot2 |> 
  filter(rowid >= 27)

#make into a list
oregon1_boot_post_list <- as.list(oregon_boot1_post[,3:ncol(oregon_boot1_post)])

oregon2_boot_post_list  <-  as.list(oregon_boot2_post[,3:ncol(oregon_boot2_post)])


#make this into a list
oregon1_boot_post_list  <- purrr::map(oregon1_boot_post_list, function(bind) (cbind.data.frame(oregon_boot1_post$rowid, oregon_boot1_post$state, bind))) |> 
  #rename the columns
  purrr::map(dplyr::rename, rowid = 1, state = 2) |> 
  purrr::map(~tidyr::pivot_wider(.x, names_from = "state", 
                                 values_from = "bind"))|> 
  #reorder to make it easier for glmnet - Oregon in the second column
  purrr::map(relocate, Oregon, .after = rowid) |> 
  #drop the rowid
  purrr::map(dplyr::select, -rowid, -Oregon)|> 
  #convert to a matrix for LASSO
  purrr::map(as.matrix)



#Do the second covariate
oregon2_boot_post_list <- purrr::map(oregon2_boot_post_list, function(bind) (cbind.data.frame(oregon_boot2_post$rowid, bind))) |> 
  #rename the columns
  purrr::map(dplyr::rename, rowid = 1,fentanyl = 2)|> 
  #drop the state
  purrr::map(dplyr::select, -rowid)|> 
  #convert to a matrix for LASSO
  purrr::map(as.matrix)


#bind both lists together
oregon_boot_post_list <- purrr::map2(oregon1_boot_post_list, oregon2_boot_post_list, cbind)


#create the post period
Oregon_post_observed <- glmnet_data |> 
  select(rowid, State, od_death_rate_half) |> 
  filter(State == 'Oregon') |> 
  filter(rowid >= 27)

Plotting ATTs

Plotting the average treatment effects on the treated (ATTs) confirms the earlier results. The treatment effect is positive for the un-adjusted synthetic control, suggesting that overdose deaths were higher after the intervention. For the adjusted synthetic control, the treatment effect is negative, suggesting that controlling for the introduction of fentanyl into the Oregon drug market indicates that the oberved overdose death rate was lower than expected given the level of fentanyl in Oregon.

#run the functions to grab the bootstrap confidence intervals
unadj_conf_int <- bootstrap_cum_effects_fixed(boot_pre_list = oregon_boot_pre_list, 
                                              boot_post_list = oregon_boot_post_list, 
                                              boot_pre_outcome_column = 2,
                                              outcome_post = Oregon_post_observed,
                                              outcome_post_column = 3,
                                              lambda_min = lambda_min,
                                              conf_int = .95)



adj_conf_int <- bootstrap_cum_effects_fixed(boot_pre_list = oregon_boot_pre_list, 
                                            boot_post_list = oregon_boot_post_list, 
                                            boot_pre_outcome_column = 2,
                                            outcome_post = Oregon_post_observed,
                                            outcome_post_column = 3,
                                            lambda_min = lambda_min_adj,
                                            conf_int = .95,
                                            penalty = penalty_factor)



#grab these and bind them to the pred table
adj_conf_eff <- adj_conf_int$effects_intervals |> 
  select(rowid, eff_low, eff_high) |> 
  rename(eff_low_adj = eff_low, eff_high_adj = eff_high)


conf_eff <- unadj_conf_int$effects_intervals |> 
  select(rowid, eff_low, eff_high) 



join_list <- list(pred_table, adj_conf_eff, conf_eff)


pred_table <- reduce(join_list, left_join ,by = "rowid")

colors_2 <- c("ATT" = "blue", "Confidence Intervals" = "lightblue")


#plot unadjusted att
plot_att_unadjust <- pred_table |> 
  filter(rowid > 16) |> 
  ggplot(aes(x = Year_half))+
  geom_line(aes(y = att, color = "ATT", group = 1)) + geom_point(aes(y = att), color = "blue") +
  geom_line(aes(y = eff_low, color = "Confidence Intervals", group = 1)) +
  geom_line(aes(y = eff_high, color = "Confidence Intervals", group = 1)) +
  # scale_y_continuous(breaks = seq(from = -.09, to = .03, by = .01),limits = c(-.09, .03)) +
  geom_ribbon(aes(ymin= eff_low, ymax= eff_high), 
              alpha=0.1, fill = "lightblue", group = 1) +
  geom_vline(xintercept= "2021_1" ,lwd=1,colour="black", linetype = 2) +
  geom_hline(yintercept = 0, lwd=1,colour="black", linetype = 2) +
  labs(title = "Synthetic Control ATT: Un-adjusted (Bootstrap CIs)", 
       y = "Att estimate",
       x = "Period",
       color = "Legend") +
  scale_color_manual(values = colors_2) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

plot_att_unadjust

#plot adjusted att
plot_att_adjust <- pred_table |> 
  filter(rowid > 16) |> 
  ggplot(aes(x = Year_half))+
  geom_line(aes(y = att_adj, color = "ATT", group = 1)) + geom_point(aes(y = att_adj), color = "blue") +
  geom_line(aes(y = eff_low_adj, color = "Confidence Intervals", group = 1)) +
  geom_line(aes(y = eff_high_adj, color = "Confidence Intervals", group = 1)) +
  # scale_y_continuous(breaks = seq(from = -.09, to = .03, by = .01),limits = c(-.09, .03)) +
  geom_ribbon(aes(ymin= eff_low_adj, ymax= eff_high_adj), 
              alpha=0.1, fill = "lightblue", group = 1) +
  geom_vline(xintercept= "2021_1" ,lwd=1,colour="black", linetype = 2) +
  geom_hline(yintercept = 0, lwd=1,colour="black", linetype = 2) +
  labs(title = "Synthetic Control ATT: Adjusted (Bootstrap CIs)", 
       y = "Att estimate",
       x = "Period",
       color = "Legend") +
  scale_color_manual(values = colors_2) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

plot_att_adjust

Key Takeaways

1). Contrary to what one might think, LASSO doesn’t always appear to select the most important covariates/predictors and analysts should be aware that a method like LASSO, which tends to almost always select synthetic control weights using the outcome variable from untreated units, can have serious limitations.

2). Analysts (including myself) probably should be more mindful of including relevant covariates in their synthetic control analysis. It is very tempting just to rely on a good pretreatment fit, but relevant covariates, especially time varying ones, can have a huge impact on treatment effects. Forcing LASSO to include causally relevant variables, even if only as a sensitivity check, can be crucial to ruling out alternative explanations.

References

Cattaneo, Matias, Yingjie Feng, Filippo Palomba, and Rocio Titiunik. 2025. “Scpi: Uncertainty Quantification for Synthetic Control Methods.” Journal of Statistical Software 113: 1–38. https://www.jstatsoft.org/article/view/v113i01 (August 2, 2025).

Pickett, Robert, Jennifer Hill, and Sarah Cowan. 2026. “Synthetic Control Misconceptions: Recommendations for Practice.” doi:10.48550/arXiv.2603.19211.

Spencer, Noah. 2023. “Does Drug Decriminalization Increase Unintentional Drug Overdose Deaths?: Early Evidence from Oregon Measure 110.” Journal of Health Economics 91: 102798. doi:10.1016/j.jhealeco.2023.102798.

Wong, Anabelle, Sarah C Kramer, Marco Piccininni, Jessica L Rohmann, Tobias Kurth, Sylvie Escolano, Ulrike Grittner, and Matthieu Domenech De Cellès. 2023. “Using LASSO Regression to Estimate the Population-Level Impact of Pneumococcal Conjugate Vaccines.” American Journal of Epidemiology 192(7): 1166–80. doi:10.1093/aje/kwad061.

Zoorob, Michael J., Ju Nyeong Park, Alex H. Kral, Barrot H. Lambdin, and Brandon del Pozo. 2024. “Drug Decriminalization, Fentanyl, and Fatal Overdoses in Oregon.” JAMA Network Open 7(9): e2431612. doi:10.1001/jamanetworkopen.2024.31612.